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Abstract 

In considering the mathematical problem of describing the geodesies on 
a torus or any other surface of revolution, there is a tremendous advantage 
in conceptual understanding that derives from taking the point of view of 
a physicist by interpreting parametrized geodesies as the paths traced out 
in time by the motion of a point in the surface, identifying the parameter 
with the time. Considering energy levels in an effective potential for the 
reduced motion then proves to be an extremely useful tool in studying 
the behavior and properties of the geodesies. The same approach can be 
easily tweaked to extend to both the nonrelativistic and relativistic Kepler 
problems. The spectrum of closed geodesies on the torus is analogous to 
the quantization of energy levels in models of atoms. 



1 Geometry as Forced Motion? 

The central idea of general relativity [TJ[2] is that motion under the influence 
of the gravitational force is viewed instead as the natural result of geodesic 
motion in a curved spacetime. It was this application to physics that gave a 
tremendous boost to interest in classical differential geometry at the turn of the 
twentieth century just as it was maturing as a mathematical discipline. The 
term 'geodesic' is not part of our common everyday language, but it is a simple 
idea that most of us have some intuition about from the role of straight lines in 
flat Euclidean geometry and great circles on a sphere: geodesies are paths which 
at least locally minimize the distance between two points in a space and are such 
that the direction of the path within the space does not change as we move along 
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it. For example, moving along a great circle on a sphere, the direction of the 
path only changes within the enveloping 3-dimensional Euclidean space to stay 
tangent to the sphere, but within the sphere itself, does not rotate to the left 
or to the right with respect to its intrinsic geometry. If one imagines an insect- 
sized toy car on a much larger sphere, moving along a great circle means locking 
the steering wheel so that the car always moves straight ahead, never changing 
direction. 

The torus (donut shaped surface) is a nice example to follow the flat plane 
and the highly symmetric curved sphere as a playground for gaining concrete 
experience with geodesies [3J, but the techniques discussed here which prove 
useful in its description apply to any surface of revolution 4-7 ], provided that 
the curve being revolved has an explicit arclength parametrization (not true of 
a simple parabola of revolution, for example) . These techniques are familiar to 
any physics undergraduate student who has taken a second course in classical 
mechanics, the motivating application which for the most part led to the de- 
velopment of calculus itself. These elementary tools are not part of the typical 
mathematician's toolbox, however, so the description of this problem one finds 
in the literature or on the web falls short of what it should be. 

One imagines tracing out a parametrized geodesic path in a surface as the 
motion of a point particle in the surface, interpreting the parameter in terms of 
which the curve is expressed as the time. Indeed a computer algebra system an- 
imation of a parametrized curve does exactly this in some appropriately chosen 
time unit. We will limit ourselves to consider surfaces of revolution about the 
z-axis in ordinary Euclidean space, so that one can describe the surface using 
an "azimuthal" angular variable measuring the angle about the z-axis (just 
the usual polar coordinate of the projection of a point into the x-y plane), and 
a radial arc length variable r which parametrizes the curve of constant angle 6 
in the surface by its arc length (not to be confused with the usual polar coordi- 
nate of the projection into that plane). Like the Kepler problem for orbits in a 
plane expressed in polar coordinates, or more generally for motion in any central 
force field in 2-dimensions, conservation of angular momentum about the axis 
of symmetry enables one to reduce the problem to that of the radial motion 
alone with the effects of the angular motion felt by an effective potential, the 
centrifugal potential, which acts as a barrier surrounding the axis of symmetry 
for motion with nonzero angular momentum. One glance at the diagram of this 
1-dimensional effective potential captures the key features of the geodesic prob- 
lem. Yet mathematicians appear to be unaware of this possibility. Visualization 
of mathematics is one of the most helpful tools for understanding, so it is worth 
bringing attention to it in this context. 

With simple modifications one can also incorporate a central force into this 
geodesic problem in 2 dimensions and apply the same machinery to reveal the 
conic section orbits of the Kepler problem, i.e., the nonrelativistic gravitational 
problem, as well as understand its relativistic generalization to motion of a test 
particle in a nonrotating black hole or outside any static spherically symmetric 
body. One does not even really need the advanced classical mechanical knowl- 
edge (the calculus of variations and Lagrangians, see Appendix B) as long as one 
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buys into the concept of the conservation of energy and of angular momentum, 
ideas already presented in high school physics. 

In this way we can reintroduce our even limited intuition about particle 
motion into the problem of finding geodesies on a surface of revolution, thus 
converting geometry back into the framework of forced motion. For an under- 
graduate who has taken separately multivariable calculus, differential equations 
and linear algebra, a little exposure to elementary differential geometry brings 
all of these tools together under one roof and has the potential to play a unifying 
role in integrating them together, something sorely lacking in many undergrad- 
uate mathematics programs for students who want or need some working knowl- 
edge of those topics. Using geodesic motion on the torus as a central example 
gives an elegant and concrete focus to this particular class of problems. As a 
classical mechanical system in physics, it is a particularly rich and beautiful 
example of motion in a central force field in 2 dimensions, a consequence of its 
rotational and discrete symmetries. 

In the present article, we suggestively introduce the main ideas needed to 
consider geodesies, but this cannot substitute for a more detailed study of the 
necessary mathematical tools. Thus it is really aimed at those who already have 
enough background, but does not exclude those with interest in the topic who 
are willing to accept the loose explanation given of the preliminary ideas. 

Finally, searching the literature one finds very little concrete discussion of 
the geodesies on the torus precisely because they cannot be fully described 
analytically, but now that computer algebra systems are readily available, one 
begins to see some limited mention of using computer algebra systems (Maple 
and Mathematica) to study them numerically. However, even if one accepts the 
defining differential equations as given and plays with numerical solutions, it 
is still important to have a good understanding of the setting in which these 
numerical games are played. The present discussion hopes to supply a detailed 
picture that is still missing elsewhere. In particular, the classification of closed 
geodesies on the torus corresponds to a discrete set of energy levels in this 
picture, mirroring the analogous quantization of energy levels in the model of 
an atom. Finding and studying these closed geodesies is equally interesting to 
those of us who enjoy mathematics for its own sake. 

2 Straight lines, circles, spheres and tori 

For straight lines in the x-y plane, the arc length of a segment is simply the 
distance between the points, which amounts to the Pythagorean theorem. To 
describe the differential of arc length between two nearby points on a curve in 
the plane, the same Pythagorean theorem is used 



ds 2 = dx 2 + dy 2 . 



(1) 



For polar coordinates in the plane 



x = r cos 9 , y = r sin 6 , 



(2) 
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since the radial and angular coordinate lines are orthogonal, one can evaluate the 
differential of arc length along a curve by applying the Pythagorean theorem to 
the two orthogonal differentials of arc length along the two coordinate lines, dr 
directly measures the differential of arc length in the radial direction, while r dO 
measures the differential of arc length in the angular direction, since increments 
of angle must be multiplied by the radius of the arc to obtain the arc length of 
a circular arc. Thus one obtains 

ds 2 = dr 2 + r 2 d0 2 , (3) 

a result which could also be obtained simply by taking the differentials of x and 
y expressed in terms of r and 9 and substituting into ds 2 , then expanding and 
simplifying the result using the fundamental trigonometric identity. 

If we repeat this in the context of Cartesian coordinates (x, y, z) in Euclidean 
space where ds 2 = dx 2 + dy 2 + dz 2 , the polar coordinates in the plane are 
upgraded to cylindrical coordinates in space 

x = r cos 9 , y = r sin , z = z , (4) 

and one passes to spherical coordinates by introducing polar coordinates in the 
r-z half plane r > 

r = psin(f>, z — pcos(j), (5) 
so that taken together one has 

x = p sin <p cos 9 , y — p sin <fi sin 9 , z^pcoscj). (6) 

Since these coordinates are also orthogonal, the iterated Pythagorean theorem 
(distance formula as a sum of squares) can be used to evaluate the differential 
of arc length along a curve, p is an arc length coordinate, while r d9 = p sin (j> d9 
continues to describe the differential of arc length along the azimuthal angular 
coordinate (9) lines, and now pd<p describes the differential of arc length along 
the polar angular coordinate ((f)) lines, so 

ds 2 = dp 2 + p 2 d<p 2 + p 2 sin 2 <t> d9 2 . (7) 

This same result could have been obtained by inserting the differentials of the 
Cartesian coordinates into ds 2 , expanding and simplifying. 
For a sphere of constant radius p — 6, this reduces to 

ds 2 = b 2 d0 2 +b 2 sin 2 <pd9 2 = [d{bcp)] 2 + [bsincp] 2 d9 2 = dr 2 + [b sin(r / 'b)} 2 d9 2 , 

(8) 

where r = b <f> is a new arc length coordinate measuring the arc length along the 
9 coordinate lines down from the North Pole of the sphere, not to be confused 
with the previous polar coordinate r. We do this to compare the sphere arc 
length formula directly with the plane formula in the same variables. If we had 
instead introduced an angle ip = n/2 — 0in the constant 9 half plane measured 
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up from the horizontal direction so that sin</> = cos ip and then a new radial 
coordinate r — b ip, we would have found instead 

ds 2 = dr 2 + [bcos(f/b)} 2 d9 2 . (9) 

Thus for both the flat plane and the sphere, a common form of the so called 
"line element" ds 2 expressed in terms of a single function R(r) which gives the 
radius of the azimuthal 9 coordinate circle corresponding to the fixed radial 
variable r within the surface 

{r , polar coords in plane, 

b^ 1 sin(r/6) , usual spherical coords (r/b = 0), (10) 
b^ 1 coa(r/b) , alternative spherical coords (r/b — n/2 — eft) 

and the standard form of the metric 

ds 2 = dr 2 + R(r) 2 d9 2 . (11) 

At this point I must remark that spherical coordinates trigger my mathemat- 
ics/physics schizophrenia. As a physicist I normally use the opposite convention 
for spherical coordinates [5] common in physics applications in which (9, eft) are 
switched with respect to the usual calculus textbook convention adopted here, 
so that (j> is the azimuthal angle in cylindrical and spherical coordinates rather 
than 9 as in polar coordinates. To make matters worse, the radial coordinates 
(r, p) are also interchanged in the physics convention! As a calculus teacher 
myself, I have to admit that it is easier to not switch the angles on students 
who have enough to keep track of as it is, and in the present context it is easiest 
to treat all of these problems with the same notation. 

To assign unique polar, cylindrical or spherical coordinates to a point in 
the plane or space, one restricts the radial variables r and p to be nonnegative 
and the zenith angle <j> to its obvious closed interval, but one has the option of 
choosing one of two obviously useful intervals for the azimuthal angle 9 



r > 
P > 

< 6 < 7T 



radial distance from z axis, 
radial distance from origin, 
zenith angle down from positive z axis, 



< 9 < 2ir or — tt < <f> < tt : azimuthal angle. (12) 

However, to describe in a continuous fashion those parametrized curves which 
wrap more than one complete revolution around the z axis, one must allow the 
azimuthal angle to take values outside these intervals. 

Both the flat plane z = within space and the spherical coordinate sphere 
of radius b arise as surfaces of revolution [IHZ] , which are constructed by taking 
a curve in the x-z plane and revolving it around the z-axis to obtain a surface 
of revolution. As long as the curve can be parametrized so that the integral 
defining the arc length parametrization can be evaluated exactly, and one can 
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invert the relation between the arc length and curve parameter to express the 
curve in terms of an arc length parametrization, we can play the same game as 
above. Revolving the x-axis itself leads to the flat plane z = 0, while revolving 
the circle x 2 + z 2 = b 2 leads to the spherical coordinate sphere of radius b. Since 
arc length is trivial for straight lines and arcs of circles, choosing any straight 
line or circle in the x-z plane will do the job. An inclined straight line leads to 
a flat cone (for example, a 4> coordinate surface in spherical coordinates, where 
r = p and R(r) — rsin^o), while a vertical straight line leads to a flat cylinder 
(for example, an r coordinate surface in cylindrical coordinates, where the new 
radial coordinate is r = z and R(r) = 1). Choosing instead a circle not centered 
on the z-axis leads to a torus. 

If the curve of this construction crosses the z-axis and its tangent line is not 
horizontal there, then the surface of revolution has a singular point where it is 
not differentiable, and the limiting tangent line if not vertical either sweeps out a 
cone with its vertex at this point, or corresponds to a limiting cusp of revolution 
if indeed it is vertical. For a torus, this conical case leads to self-intersections 
on the axis of symmetry, with two disjoint components of the surface. 

We construct a torus as a surface of revolution by revolving a circle about 
the z axis as illustrated in Fig. [TJ We take a circle of radius b > with center 
at x = a > on the x-axis, parametrized by a polar angle x measured in the 
x-z plane counterclockwise from the positive x-axis, and rotate this circle about 
the z-axis (the "symmetry axis"). We thus obtain a standard family of tori 
for which the radius of the azimuthal radius function is R = a + bcosx, while 
z = b sin x so that the Cartesian coordinates are parametrized by 

x = (a + bcosx) cos 9 , y = (a + bcosx) sin# , z = 6sinx, (13) 

or in vector form, introducing the position vector r — (x, y, z) 

r= ((a + bcosx) cos (9, (a + 6cosx) sin (9, &sinx) . (14) 

Then introducing the arc length coordinate r = bx, this becomes 

r = ((a + bcos(r/b)) cos6>, (a + 6cos(r/6)) sin#, bsm(r/b)) , (15) 

and we can add to the above list of azimuthal radius functions 

R(r) = a + bcos(r/b) , torus, (16) 

which reduces to the alternative spherical coordinate relation when a — 0. As 
before one could also derive the standard form for the line element by merely 
substituting into ds 2 the differentials of the Cartesian coordinates expressed in 
terms of the two angular variables, expanding and simplifying and rescaling the 
radial angular variable to r, but since the coordinate circles are orthogonal, it 
is just the sum of the squares of the simple circular differentials of arc length 
along the two directions. 
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Figure 1: The torus is a surface of revolution obtained by revolving about the 
z-axis a circle in the x-z plane with center on the x-axis. Illustrated here is the 
"unit ring torus" case (a, b) = (2, 1) of a unit circle which is revolved around 
the axis, with an inner equator of unit radius. The outer (\ — 0, red) and 
inner (x = ^ z7r i green) equators are shown together with the "prime meridian" 
(8 = 0, blue). The Northern (\ = 7r/2) and Southern (x = Polar Circles 

correspond to the North and South Poles on the sphere. The radial arc length 
coordinate r — b\ and the corresponding angle \ are measured upwards from 
the outer equator. The grid shown in the computer rendition of the surface 
consists of the constant 9 circles which result from the intersection of the torus 
with vertical planes through the symmetry axis (the meridians) and the constant 
r or x horizontal circles (the parallels) which result from the intersection of the 
torus with horizontal planes. 
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By introducing the shape parameter c = (a — b)/b > —1 in terms of which 
one has a/b = 1 + c and (a + b)/b = 2 + c, one finds the following distinct subsets 
of the whole family. The ring tori [5] result from the range a > b or c > 0. A 
horn torus results from a = b or c = in which the inner radius of the torus 
goes to zero, and the origin itself has a limiting behavior which is the tip of 
two cusps of revolution lacking a tangent plane so it is a bit problematic. For 
0<a<6or— l<c<0 one has spindle tori with two self-intersection points 
on the axis of symmetry and an additional part of the surface inside the outer 
surface. For the limiting case a — or c = — 1 one has the ordinary sphere of 
radius b where \ — f- Note that apart from on overall scale factor which doesn't 
change any geodesic properties of the torus (which only depend on the shape 
of the torus), we have a one-parameter family of distinct (nonsimilar) torus 
geometries parametrized by the shape parameter — 1 < c < oo, with c = — 1 
corresponding to the limiting case of the sphere. Note that the commonly used 
horn and spindle terminology is interchanged by Gray et al [5] for some reason, 
while apparently Kepler named the inner and outer surface components of the 
self-intersecting spindle tori to be lemons and apples respectively for obvious 
shape reasons [ID] . 

Note that if one expresses the differential of arc length in terms of the two 
angular coordinates and the shape parameter, one finds 

ds 2 = b 2 [d X 2 + (c + 1 + cosx) 2 rf6» 2 ] . (17) 

The overall scale factor b 2 just rescales the unit of length on the torus, but the 
so called "conformal geometry" only depends on the shape parameter c. 

To have a concrete simple ring torus example to play with, it seems natural 
to use the torus with (a, b) — (2, 1) and shape parameter c = 1, which we will 
call the "unit ring torus" (my terminology) . Any ring torus has a Northern and 
Southern Polar Circle (my terminology but obvious) like the North and South 
Poles on a sphere, and both an inner and outer equator like the equator on a 
sphere. Like the sphere, its plane cross-sections at constant 9 are circles, "the 
meridians," marking off longitude from the "prime meridian" 9 = 0, whose in- 
tersection with the outer equator is the origin of the (r, 9) coordinate system. 
Similarly each circle of revolution in the torus is a "parallel" or line of latitude 
as on the sphere. This terminology of meridians (plane cross-sections through 
the axis of symmetry) and parallels (circles of revolution about the axis) intro- 
duced for any surface of revolution also applies to the non-ring tori. The inner 
equator shrinks to a point at the origin for the horn torus and then passes to 
an outer equator of the inner lemon surface for the spindle tori. We will refer 
to the horn torus with (a, b) = (1, 1) and c = as the "unit horn torus." The 
upper/lower half of any torus above/below the equators might be called the 
upper/lower hemitorus in analogy with the hemispheres of a sphere, but one 
can also introduce the inner/outer hemitorus inside/outside the Northern and 
Southern Polar Circles. 

One can use as convenient range for either angular coordinate 9 or \ the in- 
terval [0,27r) or (— 7r,7r] but for r we will usually use the range (— bn, &7r]. One 
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can also conveniently map the pair of angles onto the unit square [0, 1] x [0, 1] 
by measuring both angles in units of revolutions (divide each angle in the range 
[0, 27r) by 2ir), identifying opposite edges of the square. This can be useful 
when one is considering closed geodesies, where this will tell us information 
about the horizontal and vertical periods of the geodesies compared to the fun- 
damental period 2n shared by both angular coordinates. Such closed geodesies 
can then be characterized by a pair of integers (m, n), given a fixed reference 
point through which they pass: the number m of vertical oscillations during n 
azimuthal revolutions around the symmetry axis (without retracing their path) . 
Note that while all of the various coordinates are periodic functions of an affinc 
parameter along a geodesic, only for a closed geodesic does a common period 
exist so that the geodesic itself as a function from the real line to the torus 
manifold is periodic. 

We will see that the inner and outer equators and the meridians are all 
geodesies, but for the ring tori on which we will concentrate our attention, the 
only geodesic which does not pass through the outer equator is the inner equator. 
With this one exception, one can classify all geodesies on a ring torus by studying 
the initial value problem at a fixed point on the outer equator. The meridian 
passing through the initial point turns out to be a geodesic. The remaining 
geodesies through this point have an obvious categorization into two separate 
families: those which cross the inner equator and thus require an unbounded 
range of values of the radial coordinate to describe (the "unbound geodesies" ) 
and the rest, for which the range of values of the radial coordinate is bounded 
(the "bound geodesies"). There are closed geodesies of each type so one needs 
three integers [m, n;p] to completely characterize them, where p — refers to 
the bound case and p = 1 to the unbound case with one exception. The bound 
terminology is standard in the physics approach to the problem to be described 
below. The inner and outer equators are closed and correspond respectively to 
the triplets [0, 1; 1] and [0, 1; 0] if we agree to associate the bound inner equator 
geodesic with all the unbound geodesies which cross it. Each meridian is a 
closed unbound geodesic corresponding to [1,0; 1]. The case [1,0; 0] does not 
exist. For the horn torus, the inner equator degenerates to a point and there are 
only bound geodesies. For the spindle tori there are instead two disjoint families 
of bound geodesies which can be labeled by p = for those confined to the outer 
apple surface and p = 1 for those confined to the inner lemon surface. In this 
latter case one must consider two separate initial value problems to describe all 
the geodesies, one on each equator. 

Although it is pretty obvious, it is worth noting that in addition to its 
rotational symmetry, each torus is reflection symmetric about the equatorial 
plane and every vertical plane through the axis of symmetry. This will also be 
true of its geodesies. In particular for initial data on the outer equator, these 
reflection symmetries of the torus correspond to reflections across the radial 
(vertical) and azimuthal (horizontal) directions in the tangent space, so it is 
enough to consider initial angles in one quadrant of the vertical tangent plane 
to study all types of geodesies which pass through that initial point. Because 
of the rotational symmetry about the vertical axis, it does not matter where on 
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the outer equator we consider the initial value problem. 

It is important to also note that the continuous symmetry group of this 
problem combined with its discrete symmetries gives the dynamical system of 
geodesic motion on the torus a rich structure making possible a great deal of 
insight that a system without symmetry does not allow. This is ultimately 
responsible for most of the analysis of the present article, and is usually missing 
in differential geometry textbooks dealing with the classical topics of curves and 
surfaces in space. 



3 Geodesic equations 

Let (x 1 ^ 2 ) = (r, 6) denote the two parameters/coordinates on the surface of 
revolution in general and on the torus in particular. The position vector as a 
function of these variables is then of the form 

f(r,9) = {x{r,6),y{r,6),z{r,9)) = R(r) (cos 9, sin 9, 0) + Z(r){0, 0, 1) , (18) 

where Z(r) is the height function with respect to the x-y plane, while R(r) is 
the horizontal distance from the z-axis. Considering each of the variables in 
turn as a parameter for a curve with the other variable fixed leads to the two 
orthogonal tangent vectors along the r and 9 coordinate lines on the surface 

e r = ^M) = fl'( r )/ cos g iSin g,o) + Z'(r)(0,0,l) , (19) 
or 

ee = = R(r) <- sin 9, cos 0,0), (20) 

where for the torus one has 

R'(r) = - sin(r/6) , Z'(r) = cos(r/6) . (21) 

The inner products of these two tangent vectors give the so called metric com- 
ponents gij = e, • €j = r, 9) whose values are given by 

g rr = R\r) 2 + Z\rf = 1 , g 00 = R(r) 2 , g r6 = g 6r = . (22) 

The first result g rr = 1 is a consequence of the fact that by choice r is an 
arc length parameter along its coordinate lines. The vanishing components 
reflect the orthogonality of the two tangent vectors, which form a basis of the 
2-dimensional tangent space. The matrix of metric components with the indices 
(r, 9) used interchangeably with (1, 2) is diagonal because of this orthogonality. 



g rr g r e \ _ ( 1 
ger gee ) \ R{rf 



(23) 



The differential of the position vector is then 

df(r, 9) =e r dr + eg d9 (24) 
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and its self-dot product defines the so called line element representing the square 
of the differential of arc length on the surface 

ds 2 = df{r, 9) ■ df(r, 9) , (25) 

analogous to the expressions in polar or spherical coordinates in the flat plane 
or in flat space. This line element is a quadratic form in the differentials of these 
parameters/coordinates whose symmetric matrix of coefficients is referred to as 
the components of the metric. Using the Einstein summation convention that 
repeated indices are summed over, and using both 1,2 and r, 9 for the indices 
as convenient, one therefore has a line element and metric component 
matrix related by 

2 2 

ds 2 = dr 2 + R(rfd9 2 = ^Y^ 9ij^^ = g^dx 3 (26) 

i=l j=l 

= g rr dr 2 + g r edrd9 + gg r d9dr + gggd9 2 , 

where 

g rr = 1 , g rB = ger = , g gg = R(r) 2 . (27) 

The inverse of the matrix of metric components has components denoted by 
g 13 , which for this diagonal matrix is just the diagonal matrix of reciprocals of 
the original diagonal entries 

S"' = l, g re = g er = 0, g ee = R(r)- 2 , (28) 

namely 

( 9 9 ee ) = {o R(r)- 2 ) ' (29) 

It is useful to normalize the two tangent vectors along the coordinate lines 

e- = 5^ 1/2 ei : e? = e r , eg = i?(r) _1 e e (30) 
to obtain an orthonormal basis of the tangent plane to the surface 

' Q<p - — - 1 — G q ' & q ^ Gf 1 ' G q — . ^31^ 

Then any vector X tangent to the torus can be expressed in the form 

X = X r e r + X e e e = X f e r + X § e § , (32) 

where the two components X r and X 6 are functions of the two coordinates. We 
can also evaluate a unit normal vector field to the surface which is the outward 
unit normal on the torus by taking the cross-product of these two orthonormal 
vectors reversed in sign 

h = -erxe § = -R(r)- 1 e r x e e = Z'(r)(cos9,sm9,0) - R'(r)(0,0,l) , (33) 
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where the radial arc length condition again confirms that this is clearly a unit 
vector. 

Since the two basis vectors e r and eg not only change in directions tangent 
to the surface as one moves around in the surface, but also must change in 
the normal direction to remain tangent to the surface as the orientation of the 
tangent plane changes, their partial derivatives with respect to r and 9 will 
consist of a part tangent to the surface and a part along the normal vector to 
the surface. Thus these derivatives can be expressed as linear combinations of 
e r and eg plus a multiple of the unit outward normal h. 

Be 

-±=T k ij e k + K ij h. (34) 

The expansion coefficients T k ij of the part tangent to the surface are called the 
Christoffel symbols. Ignoring the extra term along the normal direction, appro- 
priate if we are only interested in how vectors change within the surface, leads 
to the so called covariant derivative within the surface. This is useful since we 
want to define a geodesic curve by the property that its tangent vector does not 
change its direction within the surface (or its length if properly parametrized) 
as we move along the curve. The additional term along the unit normal is a con- 
sequence of the bending of the surface, and its coefficients arc the components 
of the "extrinsic curvature" of the surface, also called the second fundamental 
form modulo a sign convention. Here we are only interested in the intrinsic 
geometry of the surface, so we will not worry about this object, interesting in 
its own right. 

We can define the covariant derivative within the surface to be the ordinary 
partial derivative on scalar functions in the surface (functions of x l , i.e., of (r, 9), 
functions which can arise by re-expressing functions of the original Cartesian 
coordinates on space, or which are simply new functions of the surface coordi- 
nates) 

V,/ = g, (35) 

but for vector fields tangent to the surface we can simply ignore the par- 
tial derivative component along the normal direction by defining the covariant 
derivatives of the two basis vector fields using the Christoffel terms only 

Vie, = Y k ij e k . (36) 

We can then extend them to a linear combination of the two surface vector 
fields by linearity and the product rule to the following linear combination of 
products of scalars and basis vector fields 

X = X r e r + X e e e = X't* , (37) 

gvi dX* 
VjX = (VjX*)* + X\V^i) = —e t + XT fe ^e fe = ^— a + XT^a 

BX i \ 

^— e t + r jk X k ) e t = [V ej X l ] ei , (38) 
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after a convenient relabeling of the indices which are summed over. The extra 
correction terms just take into account that not only are the components of the 
vector X changing, but also the basis vectors themselves. Furthermore, this 
only measures changes of the vector X within the surface, ignoring how it must 
change in the normal direction to remain tangent to the surface. 

Since the basis vectors are orthogonal to themselves and h, we can easily 
project out the Christoffcl symbol components by taking dot products of their 
defining relation 

de- 2 
ee ■ -^j = F k ije e ■ e k = g tk T k ij = 9ekT k ij = 9ejX l % 3 (no sum on €) , (39) 

fe=i 

where we remind the reader that the repeated index k is summed over, but each 
sum only consists of one term since get = if k ^ t. It turns out that the 
Christoffel symbols are also given directly by derivatives of the metric compo- 
nents 

- 1 „u ( d 9]i , 9g k i dg jk \ 
Using either approach to evaluate them in general and for the torus, one finds 

T r gg = -i(i?( r ) 2 )' = -R{r)R\r) = (a + 6cosx)sin X , (41) 
r g r e R'ir) sin X , . 

1 r6 = 1 Br = = : T , (42) 

R(r) a + bcosx 

where we use % and r/b interchangeably as convenient. 

Suppose (x l (Xj) = (r(A), 0(X)) is a parametrized curve in the surface, leading 
to the parametrized position vector 

f(r(\),9(X)) = R(r(X))(cos 6(X), sin 0(A)) + Z(r(X))(0, 0, 1) . (43) 

The tangent vector to this curve can be evaluated using the chain rule 

-r(r(A), 6(X)) = — — + ^ = —e r + —e B , (44) 

or in an abbreviated notation 

df dr d9 dx l 

Tx = Tx er + Tx ee = ^x e - (45) 

In order to compute the derivative of this tangent vector along the curve within 
the surface, we need the product and chain rules. The derivative of the basis 
vectors along the curve, ignoring the contribution along the normal, i.e., using 
the covariant derivative, is the chain rule application 



dei dei dx J Dei y k ^ 
~dX = dx~ ~dX ^ ~dX = jl ~dX 



r fe ,i^e fc (46) 
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and defines the covariant derivative of the basis tangent vectors along the curve 
itself. The covariant derivative of the tangent vector along the curve is then a 
product/sum rule application to the following linear combination of products 

D 2 x i D ( dx l \ _ d 2 x l dx l De, _ d 2 x l dx l , dx j 

~dX*~ ei = ~d~X \~d\ ei ) ~ ~~<D? ei + ~dX~dX~ ~ ~d>? e% + ~dX jl ~d\ Ek 
d 2 x i dx j dx k \ 

ly +T ^Hx) e » (47) 

after a convenient relabeling of the indices which are summed over. 

A geodesic parametrized by an affine parameter by definition satisfies the 
following system of ordinary differential equations 

d 2 x % i dx j dx k d 2 x l v-^ ,■ dx' J dx k , 

^ + r ^A^A=°' ~dX^ + jk ~dX~dX = ' ™ 

3 = 1 k=l 

for i = 1, 2 recalling that we use k = 1,2 interchangeably with k — r,6. 
These equations are simply a statement that the tangent vector to the curve 
does not change at all within the surface as one moves along the curve. 

One can easily see that any linear change of parametrization A = a\+b where 
a and b are constants will preserve this form of the differential equations. It turns 
out that this one-parameter family of parametrizations of a given geodesic is the 
only freedom left in the choice of parameter in order that the tangent vector 
have constant length 

df df ds 2 ( ds\ 2 dx l dx^ 



dx dx-dx 2 = {dx) = g -^x^x= const - (49) 

In particular one can always introduce an arc length parametrization along the 
geodesic by choosing ds/dX = 1. Thus an affine parametrization is simply any 
parametrization for which the parameter is a linear function of some arc length 
parametrization A = As + B, where A, B are constants. 

In the present case these geodesic equations are explicitly 



r 



d8\ 2 d 2 r s „, N { d6 



«• ■ r "UJ = d>?- R ' {rm {Tx) =°- (50 > 

d 2 6 g drdfi _ (f9_ R'(r) dr d6 
dX^ + r0 dXdX ~ ~d~X 2+ R(r) dXdX 

" ^-^("M's)- - < 51 > 

where the final form of the second equation is easily verified by expanding the 
derivative. One needs to specify the initial position and initial tangent vector 
in order to determine a unique geodesic through a given point on the surface, 
aimed in a particular direction tangent to the surface. This is the geodesic initial 
value problem, appropriate for this coupled system of two second order nonlinear 
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differential equations. If one chooses A = for the initial parameter value, the 
length of the initial tangent vector determines the scaling of A with respect to 
arc length along the resulting geodesic. It is also natural to consider geodesies 
between two fixed points on the surface, the two point geodesic boundary value 
problem, but this usually has nonunique solutions, if not an infinite number, and 
numerical solution requires much more sophisticated approximation schemes 
compared to the initial value problem. Finding a shortest such geodesic is also 
a nontrivial matter. 

We can study these geodesic equations analytically (when possible) and nu- 
merically (when not). However, one can also introduce some additional quan- 
tities which enable one to more easily interpret how solutions of this system of 
differential equations behave. For those who are already familiar with differen- 
tial geometry, this can be the starting point for the remaining discussion. 

Before moving on we can note some simple solutions to these geodesic equa- 
tions without specifying R(r). For example, the 9 = 9 coordinate circles (the 
meridians on the torus) have d6/dX = = d 2 9/dX 2 , which satisfies the angular 
differential equation and reduces the radial equation to d 2 r/dX 2 — 0, with so- 
lution r = c\X + c 2 , so that r is an affine parameter on these special geodesies 
(recall it is actually an arc length parameter on its coordinate lines). Similarly 
for those points r = ro for which R'(ro) = 0, the radial equation is automati- 
cally satisfied while 9 becomes an affine parameter on these circles of revolution. 
For the torus, = R'(r) oc sin(x) leads to the values \ = 0, ±7r describing the 
outer and inner equators, which are therefore geodesies. For the ordinary sphere 
centered at the origin, the special geodesies of these two types correspond to the 
lines of longitude and the equator. For a horizontal plane seen as a horizontal 
line through the z-axis revolved around that axis, we only get the straight lines 
through its origin (which is the intersection of the plane with the z-axis) as 
these special geodesies, namely meridians. 



4 The physics approach 

So far this is the typical mathematical approach. However, let's adopt the point 
of view of a physicist, imagining tracing out a geodesic by identifying the affine 
parameter A with the time, so that the mental picture is now of a point particle 
moving on the surface, tracing out a path called the orbit of the particle, which 
refers to the unparametrized curve. The tangent vector of the geodesic curve is 
now called the velocity 

v = ^-=v r e r + v e e e (52) 
ciX 

and we might as well adopt a 2-component vector notation for components with 
respect to the basis vectors ej 

l dx % e I dr d9\ 

v= lx : {vy) = \TyTx) (53) 
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and its magnitude v = (v-v) 1 / 2 is the speed, which is just the time rate of change 
of the arc length along the curve (time rate of change of distance traveled) 

dx % dxi \ ds 



9ij lxlx) =d\- (54) 

Note that v r — v r is just the radial velocity while v e is the azimuthal angular 
velocity and v e — R(r)v is the azimuthal component of the velocity vector. 
The velocity can be represented in terms of polar coordinates in the tangent 
plane to make explicit its magnitude and inclination angle with respect to the 
radial direction within the surface. To introduce these, we simply represent the 
orthonormal components in terms of the usual polar coordinate variables in this 
velocity plane in which v r is along the first axis and v e is along the second axis 

(v f ,v 8 ) = w(cos/3,sin/3) . (55) 

The speed plays the role of the radial variable in this velocity plane, while the 
angle /? gives the direction of the velocity with respect to the direction e? in the 
counterclockwise sense in this plane. For affinely parametrized geodesies, the 
speed is constant along the geodesic; for arclength parametrized geodesies, the 
speed is 1. 

The second order differential equations for the geodesies 

lx =•••' lx =••• (56) 

may be interpreted as defining the radial and azimuthal angular accelerations 
necessary to keep the moving point on a path which doesn't change direction 
within the surface. As a system of second order differential equations, appro- 
priate initial data consists of the initial values of the unknowns and their first 
derivatives, i.e., the initial position and initial velocity. For a ring torus, every 
geodesic except the inner equator must pass through the outer equator (as we 
will see below), and because of the rotational symmetry about the vertical axis, 
we might as well fix the initial position to lie on the outer equator r = at 
= 0. It then remains to specify the initial velocity at that point, but if we 
assume an arc length paramctrization, then the velocity vector must be a unit 
vector, so only its direction within the vertical tangent plane at that initial po- 
sition remains to be specified. Thus we have a 1-parameter family of arc length 
parametrized geodesies which interpolate between the outer equator (horizontal 
initial direction) and the meridian (vertical initial direction). A computer al- 
gebra system can be easily programmed to provide a procedure for numerically 
solving the differential equations with these initial conditions as a function of 
the initial direction specified by the initial value /3 of the angle /?, which at 
this initial position is the direction measured from the vertical direction in the 
clockwise direction as seen from outside the torus looking towards the surface. 

Using such a computer algebra system program one can visualize the geodesies 
which result from various angles (3 - Decreasing the angle from /3 = 7r/2 to- 
wards zero, one sees that the path rises higher and higher before falling back 
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Figure 2: The tangent plane x — 3 at the initial position (3,0,0) at the in- 
tersection of the prime meridian and the outer equator of the unit ring torus 
where the tangent plane is vertical. The initial direction of a geodesic (unit ve- 
locity vector v = (cos /3q, sin /3q)) is here parametrized by the angle /?o measured 
from the upward vertical axis in the clockwise direction as seen from the outside 
looking in. 

towards the outer equator until at a certain point it goes over the North Po- 
lar Circle. Decreasing the angle further, the path creeps closer and closer to 
the inner equator, and then at a certain point, the path breaks through and 
"threads the donut hole," i.e., crosses the inner equator. This behavior is easily 
quantified as we will see below. 

However, an interesting class of geodesies are the closed ones which return to 
the initial position after an integral number n > of revolutions around the z- 
axis while making m > vertical oscillations, cither having a point in common 
with the inner equator (p = 1) or not (p = 0), and thereafter retrace the 
same path over and over. While the Cartesian coordinates (x, y, z) are periodic 
functions of A for all geodesies, these three functions only have a common period 
for the closed such geodesies, i.e., the position vector r is a periodic vector 
function of A. We will see below that if a geodesic returns to the initial point, 
its velocity must have either the same initial angle or its complement, and in 
the latter case, it will return to the initial point with the same initial angle after 
twice as many azimuthal revolutions, which follows from reflection symmetry 
across the outer equator. These closed geodesies can therefore be classified by 
their period number pairs and the inner equator parameter: [m, n; 0] or [m, n; 1]. 
By trial and error one can find some of these special geodesies using a numerical 
program, making it a sort of mathematical computer game. For some reason 
this fascinated me and perhaps it might do the same for some students. A Maple 
worksheet is available to test this out at my website. However, why guess if one 
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can predict? This is possible as described below. 

5 Energy and angular momentum: the key to 
unlocking the mystery 

To understand the system of two second order geodesic equations, one can 
use a standard physics technique of partially integrating them and so reduce 
them to two first order equations by using two constants of the motion which 
arise from the two independent symmetries of the "equations of motion" : time 
translation, which results in a constant energy, and rotational symmetry, which 
results in a constant angular momentum. There is no need to understand how 
this happens — it is a topic which belongs in an advanced mechanics course. Here 
we will just manipulate our specific equations to accomplish our goals. 

Since the mass m of the point particle whose motion traces out a geodesic 
path is irrelevant in this problem, those physical quantities like energy and 
momentum which involve the mass as a proportionality factor will instead by 
replaced by the '"specific" quantities obtained by dividing out the mass. Thus 
the specific kinetic energy (kinetic energy \mv 2 divided by the mass ra) is just 

1 2 1 dx % dxi 

E = —V = — QiA 

2 2 yiJ dX dX 
1 fdr\ 2 1 , , 2 fde\ 2 

= -v 2 cos 2 (3 + -v 2 sin 2 ft , 

and vice versa 

v={2E) 1/2 . (60) 

Both the specific energy and speed are constant for an affine parametrization of 
the geodesic. 

In the physics approach the specific energy of the particle is constant because 
from the point of view of its motion in space, it is only accelerated perpendicular 
to the surface (since the component of the second derivative along the surface 
has been designed to be zero). If a force were responsible for this acceleration, 
namely the normal force which keeps the particle on the surface, since it is 
perpendicular to the velocity of the particle it would not do any work on the 
particle (recall that work done is the integral of the dot product of the force 
and velocity vectors). Thus its energy and hence specific energy E must remain 
constant (said to be "conserved" along the geodesic). Similarly the speed v = 
(2E) 1 / 2 must therefore be constant along a geodesic according to this reasoning. 

In this same physics language, the second geodesic equation tells us that the 
specific angular momentum about the axis of symmetry (defined exactly as in 
the case of circular motion around an axis with radius R(r), namely the velocity 
v e = R(r) d8/dX in the angular direction multiplied by the radius R(r) of the 
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(57) 
(58) 
(59) 



circle) 



£ = z ■ (fx v) 

d9 

= R(r)v e = R(r) 2 — = R{r)v sin /3 (61) 
d\ 

is constant along a geodesic, or "is a conserved quantity" for the motion. This 
relation can then be used to re-express the variable angular velocity d9/d\ in the 
specific energy formula in terms of the constant angular momentum to obtain 
the constant specific energy entirely expressed in terms of the radial motion (i.e., 
only involving the variables r and dr/dX) and another constant of the motion 

„ 1 fdr\ 2 1 i 2 

In other words we have reduced the problem to motion in one dimension gov- 
erned only by a first order differential equation for the radial coordinate. Before 
continuing its description in the next section, some general considerations are 
worth mentioning. For simplicity we drop the "specific" qualifier from the energy 
and angular momentum. Note also that for a ring torus where the azimuthal 
radius R{r) cannot vanish, if the angular momentum £ is zero, then d9/dX must 
vanish leading to a constant value of the azimuthal angle 9 — 9q corresponding 
to purely radial motion. This is a trivial geodesic as shown above. We can 
therefore concentrate our attention on the nonradial geodesies. 

Re-expressing the velocity in terms of the two constants of the motion yields 

(v\v § ) = ((2i?) 1 /2 C0S ^ ;( 2i ?) V2 sin/3) = ^jL^j . (63) 

The second component of this equation gives a direct relationship between the 
angle of inclination and the energy for a given angular momentum for nonradial 
motion /3 ^ 0, I ^ 

f 

R{r) sin /3 (64) 



(2E) 1 / 2 



or equivalently 



E = 5— = = — . (65) 

2i?(r) 2 sin 2 /3 2i?(0) 2 sin 2 f3 

The latter relation is useful to geometrically interpret the initial velocity data at 
r = for this problem. When E = | = \v 2 , the speed v is 1 which corresponds 
to an arc length parametrization. However, we will see soon that it is more 
convenient to fix the angular momentum than the energy. For example if we 
choose instead £ = 1, then varying the initial angle /3o varies the energy, as 
illustrated in Fig. [3] Different energy levels then describe the choice of initial 
direction relative to the r coordinate circles (meridians). 

To interpret the former relation, note that dividing the specific angular mo- 
mentum (units of length times speed) by the speed leads to the specific angular 
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Figure 3: Instead of fixing the speed/energy (circles) in the initial velocity 
plane at the outer equator, fixing the angular momentum leads to a vertical 
line corresponding to all geodesies in either counterclockwise (I > 0) or clock- 
wise (£ < 0) azimuthal motion or in purely radial motion (£ = 0). Here the 
counterclockwise case is illustrated. 



momentum in a unit speed (arc length) parametrization, so that its units be- 
come just those of length: £/{2E) 1 / 2 = R(r) sin /3. In fact this is just the vertical 
component of the unit speed angular momentum vector 

1/{2E) 1/2 = z ■ (f x v) = R(r) sin (3 . (66) 

This combination of the constants of the motion is of course also constant along 
a geodesic, so that for initial data at (r, 9) = (0, 0) one has 

R{r) sin /3 = i?(0) sin /3 • (67) 

The obvious interpretation of this constant length is its value at a radial "turning 
point" r = r( cxt ) where sin/? = ±1 (see below), if it exists, where the unit 
velocity is aligned with the azimuthal angular direction and this reduces to 
±i?(r( oxt )). The existence of this constant is a consequence of the 1-parameter 
rotational group of symmetries of the torus (or indeed any surface of revolution) 
about its axis of symmetry, and indeed such a constant of the motion arises 
when the surface is invariant under any 1-parameter group of symmetries, a 
fact easily seen in the variational approach to the geodesic equations. In the 
mathematical language, this quantity is a constant discovered by Clairaut for 
geodesic motion in a surface described in a coordinate system adapted to this 
1-parameter group of symmetries [6j[7]. However, it appears that none of the 
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textbooks on differential geometry develop the tools of transformation groups or 
symmetry groups (associated with the names of Lie and Killing) , and so cannot 
remark on the interpretation available in most elementary textbooks on general 
relativity. These surfaces of revolution are invariant under the the 1-parameter 
group of rotations about the symmetry axis, and as such are "isometries" of 
the intrinsic metric. The circular orbits of this symmetry group are the integral 
curves of a Killing vector field generator, which in this case is d/d9. Translations 
in this coordinate 9 correspond to these rotations. The Clairaut conditions on 
Clairaut coordinates just state that the coordinates are orthogonal and one is 
a Killing coordinate, i.e., such that the metric is invariant under translation 
in that coordinate, and hence the components of the metric in the coordinate 
system adapted to that single independent Killing vector field are independent 
of that coordinate. The constant of the motion is just the inner product of the 
unit tangent vector (unit velocity) along a geodesic with that Killing vector field 

(v)e = v ■ — = R(r) sin(/3) . (68) 

We can use this constant of the motion to answer some interesting questions 
about the crossing angle of a geodesic with the parallels it encounters. For 
example, if a geodesic starting at the outer equator crosses the inner equator, 
it must do so at an angle which satisfies 

sin /3(A) = -^-sin/3 = ^sin/3 = ^sin/3 . (69) 
R(bTT) a — b c 

This crossing angle /3(A) is a bigger angle than /?o for a ring torus for which c > 0. 
If a geodesic starts out at the outer equator and returns to the outer equator 
R(r(X)) = i?(0), then sin(/3(A)) = sin(/?o), so its crossing angle with respect to 
the outer equator is either a) the same (3$ or b) its supplement 7r — /3q. This cross- 
ing angle cannot equal one of the two sign-reversed angles — /?o, — (tt — A)) since 
the azimuthal angle 9 is a monotonically increasing or decreasing function along 
a geodesic for nonzero angular momentum (the sign of d9/d\ does not change), 
always revolving either clockwise or counterclockwise around its symmetry axis 
as seen from above. If the geodesic returns to the initial starting point at the 
origin of coordinates, the first case a) means the initial data at the final point is 
the same as at the starting point and so it must retrace its path and the geodesic 
is closed, while in the second case b) by symmetry with respect to reflections 
across the equatorial plane, the curve must be a geodesic which will return to 
the starting point with the original initial velocity after another time interval 
of equal length A, and is therefore also closed. However, for geodesies which 
pass through the inner equator, the radial coordinate is either monotonically 
increasing or decreasing (dr/dX always has the same sign, never vanishing), and 
so cannot cross itself at the starting point with the supplementary angle and 
hence must begin retracing its path as a closed geodesic. In fact since for such 
geodesies both r and 9 are monotonic functions of A, the velocity must always 
remain within one of the four quadrants of the tangent space for which the sign 
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combinations of (dr/dX,dd/dX) are fixed. Thus the unbound geodesies either 
never intersect themselves or they are closed and retrace their own path. The 
bound geodesies may intersect themselves, but only at supplementary angles 
with respect to the meridians. 

The so called "turning points" of the radial motion occur at extreme radial 
values r = r( ext ) , where the radial velocity v r = v cos /? is zero and the velocity 
purely azimuthal. These points correspond to cos/3 = and so sin/3 = ±1 and 
arc therefore determined as the roots of the equation 

i?(r (cx t)) = i?(0)|sin/3 |, (70) 

or explicitly for the torus, r( ext ) = ±r( max ) where 

' (max ) ({<>> + b)\ sin/3 
b = X(max) = arccos I 

= arccos((2 + c)|sin/3 | - (1 + c)) . (71) 

This can only hold for the so called bound geodesies (see below) for which a 
turning point radius i?(r( ex t)) exists and hence must be larger than or equal to 
the inner equatorial radius 

R(0) | sin fa | = fl(r (ext) ) > R(bn) o -irb < r (ext) < irb , (72) 

so that sin /3o is larger in absolute value than the ratio of the inner and outer 
equatorial radii 

, . „ , R(bn) a — b c ,„„. 

sin/3 >-^777T-— " 7 = 77— ■ 73 
R(0) a + b 2 + c 

The formula for r( ext ) can be used to aim geodesies starting at the outer 
equator to hit a certain extremal radius. For example, to asymptotically ap- 
proach the inner equator one needs the initial angle to satisfy 

I sin/3 | = = — > Po = ±/? (crit) , 7T ± /3( C rit) (74) 

a + b 2 + c 

where 



/3( C rit) = arcsin yj-^j (75) 

(19.5° for the unit ring torus) defines a critical angle, or to reach the Northern 
or Southern Polar Circles the angle must be larger and satisfy 

defining another special angle 



/?( P oiar) = arcsin ( ) (77) 
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(41.8° for the unit ring torus). All initial angles smaller than this correspond to 
geodesies which cross the polar circles, while all initial angles smaller than the 
former critical angle result in unbound geodesies which cross the inner equator. 

Note that one can simply solve the energy equation describing the radial 
motion for the radial velocity and integrate 

X - X ° = L (2E-P/R(rr)^ dr - (78) 

For the torus this integral can be done exactly in terms of elliptic functions, but 
the result cannot be inverted to express r as a function of A. However, it can 
also just be integrated numerically to provide A(r), which becomes arc length 
once multiplied by the speed. Note also that at turning points of the radial 
motion dr/dX — 0, so that the denominator of this integral vanishes. Thus 
an integral for which the turning point r is an upper limit of integration is an 
improper integral, which even if it converges can still lead to some complications 
for numerical evaluation. 



6 The reduced problem of the radial motion 

The real advantage of our physics approach is in the diagram we can create 
to describe the radial motion. The second term in the energy is the "effective 
potential energy" due to the angular motion, called the centrifugal potential. If 
we only consider nonradial motion, we can fix a convenient value of the angular 
momentum while allowing the energy to be determined by the initial angle and 
therefore have a single potential energy function U describe the radial motion 

1 fdr\ 2 _ _ 1 i 2 _ 1 e 

2\dXj + _ ' U - 2R(rf ~ 2 (a + 6cos(r/6)) 2 ' (?9) 

Since the sum of the radial kinetic energy and this potential energy function 
is the constant energy for a given geodesic, once one sets an energy level, the 
motion must be restricted to the region where the difference E — U is nonneg- 
ative, with the radial velocity dr/dX necessarily zero at points where E = U , 
which determine the extreme values f7 ex t) °f the radial variable within which 
the geodesic is confined. Fig. @] illustrates one period of the periodic potential 
for the unit ring torus (a, b) = (2, 1) with the choice I = 1, where the inner 
equator has energy 1/2 and is the unstable equilibrium solution sitting at the 
peak of the potential energy graph. Geodesies with energy less than 1/2 are 
bound with radial motion confined between the symmetric extreme values of 
the radial coordinate, analogous to the negative energy elliptic orbits of the 
Kepler problem. The geodesies with energy E > 1/2 are unbound, requiring an 
infinite interval — oo < r < oo of radial coordinate values to describe continu- 
ously, and therefore the periodic potential shown in Fig. [SJ For example, the 
energy E — 3.942 corresponds to the 7 loop unbound geodesic shown in Fig. [6] 
for the initial angle /3q — 0.119 radians (6.8°). 
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Figure 4: The unit angular momentum potential (i.e., for I 2 = 1) for "non- 
radial" torus geodesies starting on the outer equator on the unit ring torus 
(a,b) = (2,1), where r = x- If the energy is less than 1/2, the geodesies do 
not make it through the hole: |r max | < n. Of course the potential is periodic 
as well, and we have only shown one period here: — tt < r < tt. For energy 
larger than 1/2, the variable r is unbounded. Five energy levels are shown: 2 
typical bound geodesies (E < 1/2), 1 unique bound geodesic asymptotically 
approaching the inner equator (E = 1/2), and 2 typical unbound geodesies 
(E > 1/2). The bound geodesies have symmetric extremal radii. The potential 
minimum corresponds to the outer equator, and the potential maximum to the 
inner equator. 

If we return to the original second order equation for the radial variable and 
replace the angular velocity using conservation of angular momentum we get 

d 2 r t f \ / \ fde\ 2 / / \ / \ ( £ V l 2 R'(r) dU , , 

This shows that the effective potential really does generate the radial component 
of the effective (specific) force through its derivative. However, we are spared 
having to deal with this second order equation since we have its first integral 
already, namely the energy relation. Note that this effective force is really 
just the radial coordinate acceleration required to keep the geodesic direction 
unchanged and is not a true force. The only real force on the motion (if this were 
a real point particle in motion) is the force normal to the surface required to 
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Figure 5: The energy levels E > | for the unbound geodesies on the unit ring 
torus of type [m, 1; 1] for m = 1, . . . , 7. Here the periodic potential merely slows 
down or speeds up their monotonic unbounded radial motion. 

keep the curve from leaving the surface, a force which arises from the constraint 
that the motion remain within the surface. For completeness we show that 
the constancy of energy implies the second order radial equation by a simple 
derivative calculation. If the energy is constant but r is not constant, then 

~ 2\dXJ + 2R^' dt ~ dX\dX 2 ~ R(rf ) " ( ' 

the expression in parentheses on the right hand side of the last equation must 
be zero, which is the second order radial equation. 

When the potential has a critical point where the derivative is zero, we 
obtain an equilibrium solution for which the radial variable is constant, namely 
a 9 coordinate line (a parallel) which is a horizontal geodesic circle of revolution. 
If the second derivative is positive (a local minimum), it is a stable equilibrium 
since for slightly larger energies the radial motion will be confined to a narrow 
interval of radii around it, while if the second derivative is negative, it will 
correspond to an unstable equilibrium since the radial motion will go far from the 
initial value. The torus has a periodic potential (we must consider r unbounded 
to describe motion which wraps multiple times around the circular cross-section) 
in the shape of a repeated "potential well" as shown in Fig. [5] for a ring torus. 
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Figure 6: An closed geodesic [7, 1; 1] for the unit ring torus with 7 loops crossing 
the inner equator, corresponding to an initial angle of f3 sa 0.119 (« 6.81°) and 
relatively high energy E — 3.942 on the unit ring torus. 



The bottom of the well where 

^ p ( » = 4 = %Tf (82) 

is at r = = \ at the outer equator (a stable equilibrium) , while the rim of the 
well is at r — ±?tb or \ — ±7r at the inner equator (an unstable equilibrium) 
where 

The "walls" of the potential act as a barrier to reaching the inner equator unless 
the energy is sufficiently high to overcome this barrier. The turning points 
r = f( max ) , 77 m i n ) of the radial motion for this case occur at the minimum radius 
R(r) allowed by conservation of angular momentum: where the energy level 
intersects the potential graph. The maximum angular velocity d9/d\ occurs at 
these turning points where the radial velocity is momentarily zero, simply by 
conservation of angular momentum since these turning points are closest to the 
symmetry axis about which the geodesic is revolving. For the horn torus, the 
rim of the potential well moves up to positive infinity, so that the nonradial 
unbound geodesies disappear. 

By fixing the angular momentum, one can consider a single potential function 
instead of a family at the cost of allowing the variable energy to describe different 
geodesies. For example, if we consider an initial point at the origin of coordinates 
on the outer equator (see Fig. [3]), we only have to vary the angle of the direction 
of the velocity to describe the 1-parameter family of distinct geodesic paths 
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which emanate from the initial position, and the different inclination angles 
correspond to the different values of the energy (equivalently of the speed). 
There are a number of obvious choices to fix I ^ in the potential function 
which only depend on the shape parameter 

p2 p2 

U = — = — 

2R(r) 2 2(a + bcos(r/b)) 2 2(c + 1 + cos(r/6)) 2 v ; 

For example, setting \£\ = R(0), R(bn),b or equivalently \£\/b = (a + b)/b, (a — 
b)/b, 1 = c + 1, c, 1 leads to values at the outer and inner equators respectively 
of 

, s 1 1 ( c \ 2 I , s 1 /c + 2\ 2 f f 

^=2' 2 {7+2) >WTW' U{bn) = 2{—) '2' 2?" (85) 

The middle choice won't work for horn tori where c = 0. For the unit ring torus 
(c = f) either of the last two choices are the same as simply setting \£\ = 1, 
which will be assumed below in the unit ring torus examples. 

This latter potential is illustrated in Figs. [4] and [5] for the unit ring torus. 
Those orbits for which |x(max)| < 7T are referred to as the bound orbits since 
the range of the radial variable is confined to a finite interval of values which 
is symmetric about the center r = because of the reflection symmetry of 
the potential about its center. In the physics language, the bound geodesies 
are "trapped" in the potential "well," while the unbound geodesies are free 
to pass from one well to the next as they continue to loop around the torus 
forever. For the energy level h, there are two geodesies. One is the inner 
equator which is an unstable equilibrium, and the other is a geodesic which 
asymptotically approaches the inner equator both from above and from below 
at the two limiting turning points of the radial motion, revolving an infinite 
number of times around the symmetry axis in each direction. If the energy 
is larger than i the orbits are referred to as unbound since one requires an 
unbounded range of the radial variable to describe them continuously. Note 
that only the inner equator geodesic does not pass through the outer equator 
(obviously!), so initial data at a point on the outer equator is enough to describe 
all other geodesies. 

Thus from the potential diagram typical of all ring tori, one sees five classes 
of geodesies: the inner and outer equators, the geodesic asymptotic to the inner 
equator, the remaining bound geodesies, and the unbound geodesies, to which 
we must add the sixth type, namely the radial geodesies since the radial coor- 
dinate lines are geodesies for any surface of revolution, as they must be because 
of the reflection symmetry about the plane cross-section of the surface whose 
intersection is the given radial coordinate line. 

On the other hand when one passes from ring tori to the horn torus by 
decreasing the shape parameter c to zero, the energy level for the inner equator 
grows to positive infinity and the nonradial unbound geodesies disappear. All 
geodesies are trapped in the infinite potential well whose walls go to infinity 
at the inner equator where R{r) = 0. Continuing into the spindle tori, one 
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Figure 7: The potential for the (a,b,c) = (1, 1,0) unit horn torus (left) on the 
interval — ir < r < tt and for the (a,b,c) — (1/2,1,-1/2) spindle torus (right) 
on the interval — 2n < r < 2?r. The horn torus has a single equilibrium solution 
at the outer equator at the origin, with infinite potential walls surrounding the 
inner equator. The spindle torus has equilibrium solutions at both the outer 
(apple) equator (center well) and the inner (lemon) equator (outer wells), and 
nonradial geodesies are confined to either the lemon or spindle. Turning points 
for one energy level are shown for both cases. 

has two separate sets of nonradial bound geodesies confined either to the apple 
or the lemon component of the torus, and the two points of self-intersection 
are protected by the infinite potential walls surrounding the axis of symmetry. 
The infinite values of the potential occur at ir^/fr = ±Xco where i?(r oc ) = 0, 
namely 

Xoc = arccos {—a/b) = arccos(— (c + 1)) . (86) 

There are two stable equilibrium solutions, one at the inner equator on the lemon 
and one at the outer equator on the apple. These potentials are illustrated in 
Fig. where the condition \£\/b — 1 is assumed so that [7(0) = l/2(c + 2)~ 2 < 
1/2 and U{bir) = l/2c -2 > 1/2. As one decreases c from towards —1 the inner 
well narrows and rises as Xoo — > n/2 and the outer well widens and deepens, 
until at the limiting case c = —1 of the sphere, the inner and outer wells are 
the same shape, and only the inner well is needed to describe the nonradial 
geodesies, all of which are bound and closed with — tt/2 < r < tt/2. Note that 
if we allow — b < a < or — 2 < c < — 1, then we again get the spindle tori 
but with r = now describing the inner equator on the lemon surface. This 
allows one to describe the lemon geodesies using initial data at the origin of 
coordinates (r,9) = (0,0) as in the remaining cases. 

Since the outer equator is a stable equilibrium, nearby geodesies will oscillate 
about it with a frequency that is easily calculated by approximating the bottom 
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of the potential well by its quadratic Taylor polynomial. The equations of 
motion for these limiting solutions are then 



d9 

0: Tx 



r= o {a + bf^ 9 (a + b) 2 ^ 



„ 1 fdrV £ 2 £ 2 r 2 
r^b: E=-\—\ + - —7 + ■ 



2 \dXj (a + b) 2 2b(a + b) 3 ' 

having imposed the initial condition 9 = at A = on the angular solution. 
Defining 

1 TZ 2 =E-— (87) 



= +-V. (88, 



[6(a + 6) 3 ]V2 ' 2 (a + fe) 2 

the energy equation describes a simple harmonic oscillator 

\ 2 

Integrating this for A as a function of r (solve for dX and integrate), imposing 
the initial condition A = at r = yields the result 

r = ±Ksm(ujX) = ±K sin((c + 2) 1/2 9) . (89) 

Thus there are (c + 2) 1 / 2 radial oscillations during one revolution about the 
symmetry axis, or y/3 rts 1.73 for the unit ring torus. The final equation here, 
having eliminated the parameter A, directly gives the orbit equation for this case, 
namely the path traced out by these limiting geodesies near the outer equator. 
Note that the smallest positive integer value of the positive shape parameter 
that leads to an asymptotically closed orbit for the small oscillations at the outer 
equator is c = 2 (i.e., a/b — 3), with 2 oscillations per revolution. For a sphere 
(limiting value c = — 1) this leads to one oscillation per revolution, corresponding 
to a great circle, which also characterizes the exact oscillations about the equator 
of any amplitude. For the redundant parameter range — 2 < c < — 1 for spindle 
tori (duplicating the range — 1 < c < 0) where r = instead corresponds to the 
inner equator (on the lemon) of the self-intersecting torus, one finds an infinite 
number of values of c for which the periodic limiting small amplitude oscillations 
have a period which is an integer multiple of the length of the inner equator. 



7 The orbits 

If we are only interested in the paths traced out by the geodesies (the orbits) , we 
can easily eliminate the affine parameter A and get a direct integral relationship 
between 9 and r in general, first solving the energy equation for the radial 
velocity and then using the chain rule 

d9 _d9 dr _ I 1 t 1 

~ dX'dX ~ R(r) 2 (2E-£ 2 /R(r) 2 y/ 2 ~ R(r) (2ER(r) 2 - P) 1 / 2 ' (90) 
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Figure 8: Periodic orbits [3, 2; 1] (left, unbound) and [3, 2; 0] (right, bound) for 
the unit ring torus corresponding to an initial angles of Pq — 0.0.3226 (18.5°) 
and /3 = 0.7167 (41.1°). 



Note that this goes infinite at a turning point of the radial motion where dr/dX = 
0. This expression can be integrated and expressed in terms of the initial data 
at (r, 9) = (0, 0) parametrized by the initial polar angle Po by using the relation 

t 2 

E= =— , (91) 

2i?(0) 2 sin 2 A, 



which leads to 

i?(0) sin f3 



o R(r)(R(r) 2 - i?(0) 2 sin 2 A,) 1 / 2 
x (c+l)sin/? 



dr = F(r, (3 ) 



o (c+ 1 + cosx)[(c + 1 + cosx) 2 - (c+ \) 2 syfl 2 fa} 1 / 2 
= G( X ,Po). (92) 

Since this integral cannot be done in closed form unless c = — 1 (the limiting 
case of a sphere), one must numerically integrate this. Introduce the reciprocal 
of the angle 9 measured in cycles (i.e., the reciprocal of 9/(2tt)) 

N(x,Po) = r( 2n R v A) < /3( C rit) • (93) 

This can be used to numerically determine the values of the initial angle /?o 
which correspond to closed orbits. 

For an unbound geodesic the value of this quantity after one radial oscillation 
is the theta frequency function N(2tt,Pq) = 2ir/G(2ir, /?o) associated with the 
theta period G(2n, Pq). This quantity is the number of times the azimuthal 
increment occurring during that complete oscillation fits into a full azimuthal 
revolution, i.e, the number of radial oscillations per azimuthal revolution. If 
this is a rational number m/n 

Tfl 

N(2n,p ) = -, (94) 
n 
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Figure 9: The number N(2tt, /3q) of radial oscillations per azimuthal revolution 
for the unbound geodesies on the unit ring torus plotted versus the initial an- 
gle fa first over the interval [0, /3( cr it)] and then for [0.9999/3( cr i t ), /3( crit )], where 
/3( cr i t ) 0.3398. The dashed line shows the asymptotic k/fio approximation for 
small angles. The rational values N(2ir,f3o) = m/n correspond to the unbound 
closed geodesies of type [m,n; 1]. 



then after m radial oscillations, n azimuthal revolutions will take place, cor- 
responding to an unbound geodesic of type [m,n;l]. Fig. |9] shows a plot of 
N(2ir,f3o) over the entire interval (3 £ [0, /3( C rit)) of values corresponding to the 
unbounded geodesies moving up initially in the positive 9 direction. The graph 
rises very quickly due to a vertical asymptote at /3q = 0, while it decreases 
extremely slowly to its limiting value limp ^p (ciit) N(2n,(3o) = where it ap- 
proaches a vertical tangent, as shown in a closeup of the endpoint behavior. 

The integral F(n, /3( cr it)) is an improper integral whose integrand denomina- 
tor goes to zero at the right endpoint, and it has an infinite value correspond- 
ing to the infinite number of azimuthal revolutions which occur as the critical 
geodesic asymptotically approaches the inner equator. The reciprocal N(tt,/3q) 
therefore goes to zero as (3q /?(crit) from the left, but very slowly. To get an 
idea of how slowly, set fto = /3( C rit) (which satisfies (c + 1) sin/3( cr ; t ) = c), and 
evaluate the Taylor expansion about \ — ir of the reciprocal of the integrand 
of F(ir, /3(crit) ) / (2tt) to find that the leading behavior of the integrand at its 
endpoint \ = 71 is 

27rc 1 /2(ir- x ) (95) 

so its antiderivative behaves like 

^H^-X)- 1 ]. 06) 
Thus if in the discrete numerical plotting of this function at the endpoint, if the 
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last sampled point is it — \ = 10 10 , say, then the value of the antiderivative 
there for the unit ring torus (c = 1) is only about 

^ln[( 7 r-x)- 1 ]«3.3. (97) 

Z7T 

The reciprocal of this gives a rough idea of the value of N reached by this last 
sampled point. This means numerically finding rational values of N for even 
relatively large proper fractions where f3g is close to the critical value will be 
difficult. For example, for the unit ring torus where /3( cr it) ~ 0.3398369094 
(19.471°), one easily finds /3 = 0.2382795502 (13.7°) for the [3, 1; 1] geodesies 
and A) = 0.0.3226432999 (18.5°) for the [3,2; 1] geodesies, but the result /3 = 
0.0.3226432999 (19.454°) for the [3,4;1] geodesic takes a considerable time, 
while the value /3 = 0.3395532232 (19.469° for the [3, 5; 1] geodesic takes too 
long for a simple root finder to obtain and one must resort to trial and error. 

The azimuthal integral 9 = G(x, /?o) f° r the orbits of the unbound geodesies 
can be easily approximated in the limit of very small initial angles where a small 
increment of the azimuthal angle occurs during one radial oscillation. Under 
the condition £/V2E — (a + b) sm(3 <C 1, the integral approaches 

bdx (98) 



f2E R{bx) 2 



which has a complicated exact value but the increment over one revolution 
— 7T < X < T is relatively simple (computer algebra system exercise) 

M= " ( ^! ft. (99) 



(a 2 - 6 2 ) 3 / 2 



This leads to 



(a 2 -b 2 ) 3 / 2 1 c 3 / 2 (c + 2 s ) 1 / 2 1 
ab(a + b) fj (c + 1) /3 

which nicely fits the approach to the vertical asymptote for small angles as 
shown in Fig. [SJ 

For a bound closed geodesic, the increment A8 in 8 during one radial oscil- 
lation cycle is four times the increment from x = to x — Xmax, so one must 
modify the number of radial oscillations per azimuthal revolution function to 

Abound) (A)) = 2lT (R Z R Z <(C+ 2) 1 ' 2 , /3 (crit) < /3 < J ■ ( 101 ) 

This upper bound for the theta frequency of the bound geodesies comes from 
the limiting case of small radial oscillations about the outer equator which is 
the largest radius circle of revolution in the torus. Discussed more at length in 
Appendix C, as the amplitude of the oscillation increases (increasing X(max))i the 
conservation of angular momentum at the smaller radii of revolution increases 
the azimuthal angular velocity compared to the radial motion, thus increasing 
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Figure 10: The number N{2tt,(3q) of radial oscillations in one azimuthal rev- 
olution versus the initial angle j3o € (/3( C rit), 7r /2], here shown for the unit ring 
torus with the upper bound 7V(2vr,7r/2) = (c + 2) 1 / 2 = V3 rj 1.73. The ratio- 
nal values N(2tt,Pq) = m/n correspond to the bound closed geodesies of type 
[to, n; 0] with to radial oscillations during n azimuthal revolutions. 

the ratio of azimuthal revolutions per radial oscillation. Fig.[T7]of Appendix B 
shows this behavior for the unit ring torus. 

The rational values Abound) (A) ) = m/n correspond to the bound closed 
geodesies of type [to, n; 0] with to radial oscillations during n azimuthal revolu- 
tions. This condition must be solved numerically for a given ratio m/n bounded 
above by the condition 

m/n < (c + 2) 1/2 . 

For the unit ring torus where (c + l) 1 / 2 = \f?> « 1.73, for example, the (in- 
equivalent) bound closed geodesies can be organized by increasing n and then 

1 1 3 1 4 5 1 3 5 1 2 3 4 6 7 8 

T ; 2' 2 ; 3' 3' 3 ; 4' 4' 4 5 5' 5' 5' 5' 5' 5' 5 ; "' ' { ' 

For any bound geodesic r is implicitly a periodic function of 9 with period 
4G(xmax(/3o), fio)- For the closed such geodesies, the above condition simply 
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Figure 11: Bound periodic orbits [1, 1; 0], [1, 2; 0] for the unit ring torus corre- 
sponding to initial angles of /3 ~ 0.4097,0.3422 (w 23.5°, 19.6°) and maximum 
radius Xmax ~ 143.6°, 173.4° and lengths 15.3,21.9. Orbits with initial an- 
gle /3q < 41.8° reach the North Polar Circle, while orbits with initial angle 
A) < 19.47° are unbound, passing the inner equator if the inequality is satisfied. 

requires that this period be a rational multiple of 2tt. 

Figs. 181 and [TT1 show the first three geodesies of this list. It is relatively easy 
to find the initial angles for the larger permitted ratios m/n for sufficiently small 
integers, but as the ratio m/n gets smaller, the limiting vertical tangent in the 
graph of iV(bound) (A>) as A) - ► Acrit) from the right creates problems for accurate 
determination of the roots where -/V^ound) is a small rational number. Because 
its reciprocal is an improper integral with a positive integrand that goes to 
infinity at the critical angle, numerical evaluation will underestimate the result 
of the integral (because numerically it will only sense some highest value of 
the integrand in the discrete approximation and miss the infinite contribution). 
Thus the true angle which gives a particular rational value to -Abound) (A)) will 
have to be a bit bigger to compensate for the lower value of the approximate 
integral. Thus one can take this underestimate as a starting angle in the full 
numerical geodesic solver and increase it slightly by trial and error until one 
closes the orbit. Even for larger values of the ratio m/n where one can accurately 
determine the corresponding initial angle, if m or n gets too large, one is of 
course limited again by numerical error in the full geodesic solver since too 
many radial oscillations or azimuthal revolutions translates into accumulated 
numerical error, so the orbit may not appear to close. When it fails to close by 
just a little, one sees a constant rate of "precession" of the node (where the orbit 
crosses the outer equator r = 0) since the crossing angle at the outer equator is 
the same each time it crosses nearby and hence it repeats the original trajectory 
but translated around the torus slightly by the amount by which the return 
position to the outer equator misses the original meridian of its departure. This 
is illustrated in Fig. [T^l 

8 Lengths of geodesies 

Since increments of the affine parameter are only proportional to the arc length 
of a geodesic segment, one must evaluate the arc length from the affine parameter 
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Figure 12: An almost closed bound geodesic close to the [7, 5; 0] closed geodesic 
(f3 = 0.6124, 35.08 degrees) with a slightly higher initial angle (0o = 0.614, 
35.18 degrees) that causes the geodesic to return to the outer equator near 
the initial meridian a bit too soon, and therefore causing the nodal point of 
return near the initial meridian to precess in the opposite direction of the orbit 
(the orbit increases the azimuthal angle, the precession decreases the azimuthal 
angle of the node). A slightly smaller initial angle would lead to precession 
in the same direction as the orbit. Trial and error adjustment of the (slightly 
smaller) approximate initial angle from the integral condition enables one to 
narrow in on the correct value to as many digits as one has patience for, subject 
to the unknown accuracy of the numerical integration scheme of course. 

increment using the constant of proportionality (2E) 1 ' 2 , which is a function 
of the energy or equivalently of the initial angle. If one has only incomplete 
information not including the affinc parameter increment, like the initial angle 
for a closed geodesic, the arc length integral remains to be evaluated, say as 
a function of the interval of radial coordinate values describing the geodesic 
segment. The corresponding affine parameter integral given in Eq. (|78p can be 
evaluated exactly in term of elliptic functions, but the resulting complicated 
formula is not very enlightening. With a slight modification, this formula may 
be used to evaluate the arc length of one complete period of the motion for a 
closed geodesic. 

To evaluate the arc length of one period of one of the unbound closed 
geodesies of the class [m, n; 1], one can first evaluate the arc length of one radial 
loop (0 < r < 2irb or < x < 2w) and then multiply by the number m of radial 
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loops. The arc length of one radial loop follows from integrating 



ds _ ds/dX _ (2E)V 2 ^ 
dr dr/dX dr/dX ' 



namely 



2,r& (2Efl 2 dr _ b f 27T Rd X 



o (2E-eyR 2 y/ 2 J Q (R 2 -e 2 /(2E)y/ 2 



2 k 



Rd X 



L = 4b — ^ -, (105) 



= b / ^5 . (104) 

Jo (i? 2 -i?(0) 2 sin 2 /3 )i/2 ^ > 

For example for the 7 loop geodesic with energy E = 3.942 on the unit ring 
torus, the length of one loop is 6.45 and of the entire 7 loops: 45.12. This 
compares to the circumferences of the inner equator, the Polar Circles and 
the outer equator: [2tt(1), 2tt(2), 2tt(3)] w [6.28,12.57,18.85]. These integrals 
are easily done numerically with a computer algebra system. Notice that for 
very high energy (very small initial angle /3q — > 0), this reduces to the radial 
circumference of the torus. 

For any bound geodesic, one can easily evaluate the length of one transit 
departing from the initial point on the outer equator, traveling to a given turning 
point, returning to cross the outer equator and going to the other turning point 
and back again by taking 4 times the length from the outer equator to one 
turning point 

f x (m ax) Rdx 

{R 2 - R(0) 2 sin 2 p ) 1/2 
where 

X(max) = arccos ((2 + c)| sin A)| - (1 + c)) . (106) 

For example the [1,1 ;0] closed geodesic shown in Fig. [TT]on the unit ring torus 
has length 15.26. 

One can also categorize the closed geodesies by the number of self-intersec- 
tions they exhibit, as initiated by Irons [3|. To examine geodesic crossings 
without trying to follow a geodesic around the 3-d torus in rotatable computer 
graphics images, it is useful to map the torus onto the unit square by introducing 
two angles measured in cycles rather than radians 

(H,9)= (|^) , withO< X ,0<27T— ►O<H J 0<1. (107) 

For comparison the grid of all positive pairs [to, n] describing the distinct 
closed geodesies of the flat torus is briefly discussed in the Appendix. Like the 
flat torus, the [m, n;l] closed geodesies on a ring torus exist for all positive 
pairs (to, n) and have no self-intersections, but the [to, n; 0] closed geodesies are 
constrained by the theta frequency condition m/n < (c + 2) 1 / 2 of the limiting 
geodesies near the outer equator, which have the largest period of the bound 
geodesies, and self-intersections are allowed. As one gets closer and closer to 
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the critical geodesic which asymptotically approaches the inner equator, the 
closeness to the unstable equilibrium there allows an arbitrarily large number 
of azimuthal revolutions near it before the geodesic flies off towards the outer 
equator again, which can be arranged to happen to hit any point on the outer 
equator. Thus one can get all pairs (m,n) for m/n 0. Numerically this is 
problematic, however, since numerical error is hard to control near an unstable 
equilibrium. 

Given that only the bound geodesies can have self-intersections, comparison 
of the low integer geodesies of type [m, n; 0] reveals the following apparent cross- 
ing rules as a function of the number m of radial oscillations. For the bound 
geodesies on the unit ring torus starting at the origin of coordinates with initial 
angle j3 £ (/3( C rit)> 7r /2) orbiting in the increasing 9 direction, looking down from 
above at the upper hemitorus one finds by inspection of many cases either 0, 1, 
or 2 equally spaced series of self-crossings at 0, 1, or 2 distinct values Xx of the 
radial angle \ € [0,X(crit)]- The equal spacing increment is 2w/m. 

n = 1: No self-crossing point radii. 

n = 3: 1 self-crossing point radius: Xx ■ The series is shifted right from 8 = 
by (27r/m)/4. 

n > 3 odd: 2 self-crossing point radii: Xxi < Xxi- The 2 series are shifted 
right /left from 9 = by (2%/m)/4. 

n = 2: 1 self-crossing point radius Xx = 0. 

n > 2 even: 2 self-crossing point radii: first series at Xxi = 0, second series 
Xxi shifted right from 9 = by (2n/m)/2. 

It would be nice to understand how these rules come about, and to evaluate the 
angles Xx at which these crossings occur. 

Finally while it is convenient to describe a geodesic as the curve between 
two points which minimizes the distance between them on the surface, it is 
also important to recognize that the converse is not true. While a geodesic 
connecting two points is a critical point of the arc length function on the space of 
paths between the two specified points, as follows from the variational approach 
to this problem, there can be many geodesies of differing lengths connecting any 
two points on a surface. In the flat plane, the geodesic is the unique straight 
line between two points, but on a sphere although there is a unique great circle 
containing any two non-antipodal points, there is still a choice of the longer or 
shorter segment of the great circle to follow from one point to the other, and 
for antipodal points there are an infinite number of great circles to choose from. 
For the torus the situation is even more extreme: through any two points on the 
torus, there are an infinite number of both bound and unbound geodesies which 
connect them. However, for most point pairs there is still a minimal length 
geodesic connecting them. While one can try to aim a geodesic departing from 
the first point to hit the target point using trial and error, it is interesting that 
completely different approaches can be taken to numerically converge on this 
connecting path [I5HT7] . 

An interesting illustration of the minimal geodesic question regards the min- 
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imal length geodesic which connects two distinct points on the outer equa- 
tor. According to the approximate analysis above, geodesies departing from 
(r, 9) = (0, 0) with an initial angle differing very little from f3 = ±ir/2 oscillate 
about the outer equator with half period A9 = 7r/(c + 2) 1 / 2 , thus returning 
to the outer equator first at 6 = A9. Examining numerical geodesies for such 
initial conditions then shows that for the actual geodesies, as the amplitude of 
the radial oscillation increases, the actual half period increases slowly, so that 
the return point migrates slowly beyond 9 = A9, but the length of this geodesic 
segment remains smaller than the length of the direct route between the end- 
points along the outer equator. Thus for two points on the outer equator with 
separation equal to or less than the half period A9, the outer equator segment 
between them is the shortest geodesic between them, but for points with separa- 
tion greater than this critical length, the shortest geodesic is instead a half radial 
oscillation connecting the two points. Antipodal points on the outer equator are 
the extreme case of this type, where the [1, 1; 0] geodesic which connects those 
points is significantly shorter than half the circumference of the outer equator. 

The horizontal "focusing" of the bound geodesies which emanate from the 
origin of coordinates on the outer equator, like the crossing of great circles on 
the sphere at antipodal points, is a consequence of the positive curvature of 
the torus on the outer ring between the Northern and Southern Polar Circles. 
On the inner ring between these Polar Circles, the curvature is negative and 
geodesies from the outer equator which reach this region instead diverge from 
one another before crossing on the inner equator due to the topology of the 
torus. These geodesies begin crossing each other at the critical length L r = nb 
equal to half the radial circumference of the torus, compared to the azimuthal 
critical length Lb = (a + b)A9 = 7r(a + b)/(c + 2) 1 / 2 , leading to the ratio 
Lg/L r = (c+ 2) 1 / 2 . Irons [5] has evaluated the various curvature quantities on 
the torus, but curvature is a long story that cannot be told here, but which is 
described in Appendix C. 

9 Kepler problem? 

The beauty of mathematics is that many seemingly different problems actually 
share the same mathematical foundation and hence can be treated with the 
same methods to understand them. From geodesies on these special surfaces 
of revolution, by relaxing the condition g rr = 1 to g rr = f(r) we can easily 
consider general surfaces of revolution for which the radial coordinate does not 
measure arc length [5] . This enables us to handle the simple parabolic surface of 
revolution z = x 2 + y 2 = r 2 , for example, where g rr = 1 + 4r 2 . However, we can 
also handle the simplest example of a surface of revolution, the flat plane z = 0, 
but in the framework of polar coordinates, where resolving the geodesic equa- 
tions to reproduce straight lines is complicated but doable. Polar coordinates 
are simply not adapted to the geodesic structure of the flat plane. The real 
payoff, however, is to consider the simple addition of an external rotationally 
symmetric force to these same flat plane equations, enabling us to handle the 
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Figure 13: A typical Kepler potential for the nonrelativistic (Newtonian grav- 
ity) and relativistic (general relativity) cases, showing the "centrifugal barrier" 
near the source (small r). Circular orbits occur at the extrema of the poten- 
tial. Bound orbits have energy E < 0, while unbound orbits have E > 0. The 
relativistic potential turns over close to the source where stronger gravitational 
attraction overcomes the centrifugal force. The lowest two energy levels shown 
correspond to bound orbits with a minimum and maximum radius, while the 
upper two are unbound, but in the relativistic case, the highest energy level 
orbit overcomes the centrifugal barrier and is captured by the source, while the 
next lower energy orbit is reflected from the centrifugal barrier at the minimum 
radius in both cases (the outer portion for the relativistic case); the relativistic 
case also contains "trapped" orbits inside the barrier. 

nonrelativistic Kepler problem with the same approach, and indeed the closely 
related relativistic Kepler problem of planar motion in a spherically symmetric 
vacuum gravitational field arising either from a spherical body or a black hole. 
In both cases, the polar coordinates are essential to resolving the equations of 
motion. 

For the flat plane case R(r) = r and the energy equation becomes 



To incorporate an actual rotationally symmetric radial force F r [r) = —dU r (r) j dr 
into the mix, we just need to add its potential U r (r) to the effective potential 




(108) 
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term already present 



i (dry ie 2 . s 

E =2{Tx) + 2^ + U ^- ^ 

The new total potential then governs the problem in the same way that the 
effective potential alone did the job in our previous discussion. With a potential 
of the form 

tt r \ h k 2 , . h 3fc 2 , 11ri , 

u r( r ) = ~— ~ J3 ~> F r( r ) = ~ -^a > ( n °) 

we can handle the nonrelativistic Kepler problem for motion in an inverse square 
force if ^2 = (attractive like gravitation if k\ > 0) or the relativistic Kepler 
problem (ki > 0, ki > 0, see dJ), see Fig. Q2J The additional inverse cubic 
potential term in the latter case corresponds to an additional inverse quartic 
force that is responsible for the precession of the orbits which in its absence are 
conic sections. 

Amusingly enough, the radial oscillations about the outer equator of the 
torus are directly mirrored in the Kepler problem, and the frequency one cal- 
culates for the oscillations about a stable circular orbit is called the epicyclic 
frequency, an important astrophysical quantity. While these more general prob- 
lems are as easily dealt with as the torus geodesies, they are perhaps best left to 
a more comprehensive discussion elsewhere |12j . However, a slight adjustment 
of the Maple worksheet for the torus allows experimentation with these kinds 
of orbits as well. With a little more work, one can include a light pressure drag 
force to model the effect of stellar radiation pressure on orbiting dust particles, 
the so called Poynting- Robinson effect [13] . 

Perhaps the most important concept highlighted here that is missing in most 
undergraduate mathematics programs is energy. There is usually just no time in 
an introductory course on differential equations to introduce it (and the standard 
textbooks do not mention it), so whatever vague idea mathematics students 
have of it from high school or possibly from college freshman physics is never 
connected concretely to differential equations and the natural way it arises there. 
Even in the present article, we have danced around this connection without 
directly confronting it cither, but hopefully its utility at least in this problem 
is clear. The usual approach to geodesies through a variational approach is 
the natural way to make this connection |14j . This is an interesting advanced 
topic either for an undergraduate mathematics or physics setting. For advanced 
undergraduate mechanics in the physics setting, this approach to constrained 
motion (motion on a surface) is a natural gateway for understanding the tools 
of general relativity, as developed by Walecka [2] , for example. 



10 What next? 

It would be sad to shine the light on an interesting problem and completely 
exhaust its potential for entertainment. There is much left to explore. On the 
mathematical side, the spectrum of closed geodesies on the torus depends on the 
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shape parameter. Passing to the horn torus c = from the ring tori, the local 
maximum in the periodic potential at the inner equator moves off to infinity, 
and the unbound nonradial geodesies disappear. All the nonradial geodesies are 
repelled from the single point inner equator by this angular momentum barrier. 
For two-point self- intersecting tori — 1 < c < 0, this infinite angular momentum 
barrier splits into two barriers surrounding the two self-intersection points on 
the symmetry axis, leading to two disjoint families of nonradial bound geodesies, 
those confined to the outer apple surface and those confined to the inner lemon 
surface. For the shape parameter c = 2 (for example: (a, b) = (3, 1) so the 
ratio of the outer equator radius to inner equator radius is 2), exactly two radial 
oscillations about the outer equator fit into one azimuthal revolution in the limit 
of zero amplitude oscillations. This means that for nonzero amplitudes where 
the exact frequency will be slightly different from 2, the nodal points where the 
oscillations intersect the outer equator will slowly migrate or "precess" in the 
physics language. What is the precession frequency? This can be analyzed by 
keeping the next term beyond the quadratic one in the Taylor expansion of the 
potential about its minimum. (Similarly one can try to analyze the precession 
of the node for near miss closed geodesies to calculate the incrementing angle 
of precession described above.) For the pinched torus, all radial geodesies pass 
through the pinching point at the origin, but all of them have a vertical tangent 
vector, so one looses uniqueness of the initial value problem for the geodesies 
there; one can match up radial geodesies of different azimuthal angle on either 
side of the pinching point. What does this do to the torsion and differentiability 
of those curves? 

On the physical side, one can easily explore motion in various classical me- 
chanical force potentials as hinted at for the Kepler problems. In an advanced 
mechanics course, where only formula manipulation was possible when I was a 
student, there are exciting possibilities for student experimentation. 

Finally the obvious must be pointed out. None of this would be possible 
without the incredible tool now at the disposal of faculty and students alike: 
the computer algebra system on our personal computers. Both of the two leading 
such software packages Maple and Mathematica offer to motivated students an 
enormous potential for exploring mathematical and physical problems, if given 
a little guidance from faculty. The problem for undergraduate teaching is to 
find the opportunities to engage students in picking up this tool in parallel with 
their gradual increase in theoretical sophistication in either a mathematics or 
physics setting. 

A The flat torus closed geodesies 

For comparison it is useful to consider the spectrum of closed geodesies for the 
simpler flat torus. The simplest realization of the flat torus is the plane with 
points identified if their Cartesian coordinates differ by integers, thus reducing 
distinct points to the unit square in the first quadrant with opposite boundary 
line segments identified. One can classify these geodesies by the pair of periods 
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Figure 14: The [2,3] closed geodesic through the origin before and after im- 
posing the periodic boundary conditions for the flat torus on the unit square. 

[to, n] of the number of revolutions in the horizontal and vertical directions 
respectively. Fig. [T4l shows a typical geodesic of type [2,3] through the origin 
returning the origin at the point (2, 3) in the plane, with length V2 2 + 3 2 = \/l3. 
Fig. [15] shows in the same way a representation of the lattice of all the distinct 
closed geodesies through the origin of type [to, n] for to, n < 7, with length 
\/m 2 + n 2 . 

One can reinterpret Fig.[T5]as a representation of the actual closed geodesies 
in the tangent space at the origin using the exponential map. This map asso- 
ciates a point in the tangent space at a point P with a point Q on a geodesic 
through it by assigning to it a tangent vector with the same direction as the 
tangent vector to the geodesic at the point P, but with a length equal to the 
arc length from P to Q along that geodesic. One can therefore reproduce this 
figure for the bound and unbound geodesies on the curved unit ring torus. One 
only needs to assemble the pairs (/?o,L(/3o)) of initial angles and corresponding 
lengths for the low integer (m, n) pair bound and unbound closed geodesies and 
plot them in the vertical tangent plane at the origin of coordinates. Fig. [TS] 
shows this result. One can read about the exponential map in most mathemat- 
ical textbooks which discuss geodesies |18) . 
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Figure 15: The closed geodesies on the flat torus on the unit square charac- 
terized by the period pairs [m, n] for m,n < 7 before imposing the periodic 
boundary conditions. 



8O-1 / 




Figure 16: The closed geodesies on the unit ring torus with initial angles 
< Pq < tt/2 characterized by the triplets [m,n;p] for m,n < 7, p = 0,1, 
represented in the tangent plane at the outer equator by a line segment of 
length equal to the arc length of the geodesic and with the same angle (3q 
measured clockwise from the positive vertical axis. Included are the vertical 
meridian and the horizontal outer equator of lengths 2n and 6n respectively. 
The unbound (top) and bound (bottom) geodesies are separated by the critical 
angle /3q = 0.3398369094 (19.5 degrees) indicated by the dashed line. 
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B Variational approach to the geodesic equa- 
tions 



One can side step all of the Riemannian geometry discussion of covariant dif- 
ferentiation and parallel transport by using the variational approach which pro- 
duces the geodesic equations by extremizing an action functional on the space 
of all curves connecting any two fixed points on the torus. One can either ex- 
tremize the arc length of the curve, which is the integral of the length of the 
tangent vector in any parametrization of the curve, namely the speed function 




Actiani = / ds = / — dX = / \ — + -7T dX > ( m ) 



which is clearly independent of a change of parametrization, or extremize the 
integral of half the length squared of the tangent vector (the factor of two is for 
convenience) 

^-U(*)'*-U(£)' + *(£)' A - (ll2) 

which is equivalent to the previous case only for affinely parametrized curves for 
which the speed ds/dX is constant. These facts require a separate study of the 
topic called the calculus of variations [TS] . The integrand is called a Lagrangian 
function, and it is a function of the curve and its tangent vector. The second 
Lagrangian function is just the energy function 

/ „ dr d6\ lfdr\ 2 1 , / d6\ 2 , s 

L2 {^dX->dx)=2{dx) + 2^U) = E (m) 

while the hrst Lagrangian L\ = ds/dX is the speed function 



M^.£W(f)W(f)'- <-» 

Both are independent of the azimuthal angle because of the rotational invariance 
of the problem 

r)T 

R = R(r)^^=0. (115) 
The Lagrangian equations are a consequence of extremizing the action 

d f OLj \ dLj 

dX\d(dx3/dX)J dap' K ' 

For the second Lagrangian this produces the geodesic equations used in the main 
body of the article, with the angular equation directly giving the constancy of 



44 



the angular momentum I = dL2/d(dO/dX). Using instead the first Lagrangian, 
one obtains the constancy of the momentum conjugate to 9 



Po 



dL x R 2 d9/dX 

iO/dX) 2 

(117) 



d{dO/dX) ^/(dr/dX) 2 + R 2 (dO/dX) 
R 2 d9/dX _ £ 



ds/dX \/2E' 

where the last line holds not only for an affine parametrization but for any 
paramctrization since the expression is independent of the parametrization, and 
the final ratio of constants is just the constant vertical component R(r) sin /3 of 
the unit speed angular momentum by Eq. (|66p . For the bound geodesies where 
| sin/3| = 1 at the radial turning point r = rr ex t), its absolute value equals the 
azimuthal radius at the extremal radius: \pg\ = -R(f( ex t))- 

We can invert the first line of Eq. (|117|) , solving it for the angular velocity 
d9/dX as a function of the canonical momentum pg and then re-express the 
Lagrangian function L\ in terms of that momentum rather than the angular 
velocity 

d9 pgdr/dX ds R\dr/dX\ , . 

tx = ± r-jw=^> i = dx = ^m=^- ( 8) 

Suppose we consider geodesic segments for which r is a monotonically increasing 
or decreasing function between r\ and T2 so dr/dX ^ 0, allowing it to be used 
as a parameter itself along the geodesic, i.e., segments not containing any radial 
turning points except possibly at the endpoints. Then one can integrate the 
previous equations with r = X 

7= dr = ±F(r 1 ,r 2 ;pe) , 
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S21 = I -j:dr= / dr = G(ri,r 2 ;pe)- (119) 

-Pi 

These integrals can be evaluated exactly in terms of elliptic functions, the first 
as a long combination of such functions, while the second is a relatively short 
expression in one such function. The first of these two gives the orbit equa- 
tion 9 = 9(r) for a given specified momentum pg. This relationship cannot be 
inverted to give r = r(9). However, it can also be thought of as a very com- 
plicated condition on the momentum that in turn determines the initial angle 
that a geodesic has to be shot from its initial point (ri,0i) to terminate at 
the final point (r^,^). Given that solution for pg one can then evaluate the 
arc length of the corresponding geodesic using the second integral formula. For 
near enough points on the torus, there will be a unique shortest length geodesic 
which satisfies the previous condition, but for most bound geodesic segments, r 
is not monotonic for the shortest length geodesic between two given points, but 
involves one or two radial turning points. 

This integral function F, say when < r± < r2 < for concreteness, is 
a monotonically increasing odd function of the momentum pg on the closed 
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Figure 17: The spray of geodesies leaving the origin of coordinates in the first 
quadrant on the unit ring torus for — it <r<Tr,0<9< 2ir. As one increases 
the momentum pg from 0, first one passes through the unbound geodesies which 
cross the inner equator at r = ir (up to the thick dashed curve which asymptot- 
ically approaches the inner equator), then the bound geodesies with a smaller 
and smaller turning point radius r( max ) . Looking at where these geodesies cross 
the Northern Polar Circle at r = tt/2, starting at 9 = the crossing point moves 
to the right through unbound geodesies and then into the bound geodesies where 
their turning point lowers until it reaches that circle (the thick black geodesic, 
very close to the [3,2,0] closed geodesic which is slightly higher), after which 
the geodesies which reach that circle overshoot it first, crossing over and then 
returning to that circle as their second crossing point moves to the right and 
the turning point rises. Note also the half wavelength A9 = tt/v3 « 1.81 
of the small oscillations about the outer equator. For < 9 < A9 only the 
outer equator itself from this spray reaches points on that outer equator, but 
for A9 < 9 < 2A9 on the outer equator, a second member of this family reaches 
the outer equator (with shorter length) , while for 2 A9 < 9 < 3A9 a third mem- 
ber of this family reaches the outer equator. Of course for ir < theta < 2ir, 
shorter geodesies arrive from the opposite azimuthal direction. 
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Figure 18: The plot of Ad — 82 — 6\ versus the canonical momentum pg 
for the unit ring torus for geodesies from r\ — to r 2 = ir/2 (outer equator 
to the Northern Polar Circle), illustrated in the previous figure with 61 = 0. 
Increasing pg from through the unbound geodesic directions to the bound 
geodesic directions until the maximal turning point radius r decreases to 7r/2, 
62 — 0\ increases to the vertical tangent at pg — 2 = R(w/2) and then the curve 
turns over as the maximal radius again increases while pg decreases to 1 = R(ir) 
where a vertical asymptote occurs corresponding to the asymptotic approach to 
the inner equator. 

interval [— i?(r 2 ), i?(r 2 )] with value at pg = and which is concave up when 
Pe > 0j having an extreme value at each endpoint where a vertical tangent occurs 
when the radial turning point decreases to r 2 as \pg\ increases to i?(r 2 ). At this 
point one can decrease \pg\ from -R(r 2 ) while allowing the radial turning point to 
occur within the geodesic segment as that segment extends farther and farther 
past the original limiting angular interval, corresponding to the curve 62 — 81 
versus pg turning back on itself after the vertical tangent. Corresponding to this 
new situation, the geodesic segment requires overshooting, i.e., going outside the 
interval from r\ to r 2 to reach a radial turning point value r( ext j beyond r 2 and 
return back to it, and the discussion gets much more complicated. See Fig. [P7l 
One must then break up the azimuthal angle integral into two separate integrals 

/"• (max) jq /-r (max) ,n 

92-9!= T dr+ T dr > ( 12 °) 

J ri (XT J T2 

where the relation i£(rr max )) = \pg\ is inverted to give r( max ) as a function of 
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be I 

R(r(max)) = a + ^cos(r (max) /6) = \p g \ , r (max) = b arccos 0— L - V (121) 

Now as one decreases from the value R(r2), r( max ) increases from r2 and 
I #2 — #i| continues increasing as shown in Fig. [T8]for Pe > 0. 

This complicates the problem by making the upper limits of integration de- 
pend explicitly on the momentum as well, so that it becomes a more complicated 
condition for that momentum. The two separate integrals correspond respec- 
tively to the rise of the radial coordinate to the radial turning point and then 
its subsequent fall to the final value of the radius. These are also both improper 
integrals since the denominator goes infinite at the radial turning points, mak- 
ing the problem more tricky numerically as well to find a numerical value of pg 
which will satisfy the condition. The worst case would involve geodesic segments 
from the upper hemitorus to the lower hemitorus which require overshooting on 
both ends of the curve. 



C Convergence of geodesies 

This is an aside for the experts. We derived the orbit equation for small geodesic 
deviations from the outer equator, which can be reparametrized by the arclength 
s = (a + b) 9 along the outer equator 

r = ±TZsm((c + 2) 1 / 2 9) = ±Hsm{w s s) . w s = ^' + ^ = b' 1 (c + 2)' 1 / 2 . 

a + b 

(122) 

These oscillations have the azimuthal angular half-period Tg = n/ (c+ 2) 1 / 2 and 
the arc length half-period L — (a + b)T e = b(c + 2) 1 / 2 n — tt/u} s . Note that 
L > &7r (the half-circumference of a meridian) and L < (a + b)ir — b(c + 2)tt 
(half-circumference of the outer equator) for all allowed values of c > — 1 except 
c = — 1 which corresponds to the limiting case of a sphere where all three of 
these lengths coincide. The limit c — > oo, for which L — > oo, corresponds instead 
to opening up the torus to an infinite flat cylinder, 

This result is just a special case of the geodesic deviation equation [20] 
evaluated for deviations from the outer equator geodesic. As illustrated in 
Fig. [TS1 we have seen that points on the outer equator separated by a distance 
greater than L can be connected by a shorter geodesic which deviates from 
the outer equator. This is an example of the focusing of geodesies by positive 
curvature. The geodesies which emanate in all directions from an initial point 
on the outer equator begin crossing each other first at this distance along the 
outer equator. 

One can evaluate the single independent nonzero orthonormal component of 
the intrinsic curvature tensor (just half the curvature scalar in 2 dimensions [21] 
and equal to the Gaussian curvature of the surface [22]) to find [3] 

f _ cos(r/b) 
U efe ~ b(a + bcos(r/b)) ■ 1 6) 
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For ring tori where a + bcos(r/b) is always positive, this is negative on the inner 
hemitorus, zero on the Polar Circles and positive on the outer hemitorus, whose 
maximum (positive) value occurs at the outer equator 



^l-o - - , (124) 

This is where geodesies cross more quickly than at any other location on the 
torus, i.e., those which remain close to the outer equator will experience the 
maximum convergence. This convergence of geodesies by regions of positive 
curvature leads to the physical application of gravitational lensing. 

The unit vector e? = e r is parallel transported along the outer equator 
(as are eg and e§) while remaining orthogonal to the outer equator. With the 
geodesic separation vector £ = r(s)ef 7 the geodesic deviation equation reduces 
to 

^=-iT^| r=0 r = -w s 2 r, (125) 

which leads to the same solutions (I122j) that we started with. Thus we can con- 
clude that for any two points on a nontrivial torus (c > —1) separated by less 
than the radial half-circumference bn which in turn is less than the convergence 
arclength L, there will be a unique geodesic of lesser arclength which connects 
them. Fixing our attention on a given point of the torus, a Euclidean sphere 
in space of radius b will therefore enclose a region of the surface on which the 
geodesic from this central point to any other point in this region will be unique. 
In other words by evaluating the Euclidean distance between two points on the 
torus, one has a sufficient condition to guarantee the uniqueness of the bound- 
ary value problem for the geodesic equations. Both leading computer algebra 
systems Maple and Mathematica have numeric ordinary differential equation 
solvers which can now handle boundary value problems, so they will produce 
a geodesic between two specified points. Thus when the two points are within 
a Euclidean sphere of radius 6, these programs will deliver the minimal length 
geodesic. 

Finally for spindle tori, the curvature scalar approaches negative infinity as 
one moves in along a meridian towards to symmetry axis where a + bcos(r/b) = 
0, and then again turns positive in the lemon part of the surface decreasing from 
positive infinity to a minimum positive value at the inner equator on the lemon. 
The curvature scalar there has the positive value 

R r @fS \ r =b* = ^ = ^z^y > ( 126 ) 

which is larger than the value at the outer equator. The associated convergence 
length 6(— c) 1 / 2 ^ is therefore shorter on the inner equator (but still longer than 
the half-circumference of the inner equator), but increases as c decreases to 
the limiting value —1 where the inner and outer equators merge into the single 
equator of the sphere. 
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This is a standard topic in any upper level undergraduate treatment of 
classical mechanics. 

[20] Geodesic deviation is used in general relativity for interpreting the effects 
of curvature on the motion of test bodies, and hence is discussed in most 
textbooks on general relativity, for example, C. Misner, K.S. Thorne and 
J. A. Wheeler, Gravitation (Freeman, Sanfrancisco, 1973). 

[21] The curvature tensor emerges as soon as one considers second covariant 
derivatives and is discussed in all textbooks on differential geometry. See: 
: / / en . wikipedia.org/wiki/ Rie mann_curvature_tensor| 

[22] The Gaussian curvature is the determinant of the component matrix of the 
mixed extrinsic curvature tensor K l j = g lk Kkj (also called the shape oper- 
ator with components S l j) and is discussed in all textbooks on differential 
geometry (for example, Gray et al [5] above). See: 
http: / /en. wikipedia.org/wiki/Gaussian_curvature 
http://mathworld.wolfram.com/ShapeOperator.html 
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